Computing intersections with polygons¶

Binder IPYNB HTML

Clipping and intersection functions can be used to extract trajectory segments that are located within an area of interest polygon.

In [1]:
import pandas as pd
import geopandas as gpd
import movingpandas as mpd
import shapely as shp
import hvplot.pandas 

from geopandas import GeoDataFrame, read_file
from shapely.geometry import Point, LineString, Polygon
from datetime import datetime, timedelta
from holoviews import opts

import warnings
warnings.filterwarnings('ignore')

opts.defaults(opts.Overlay(active_tools=['wheel_zoom'], frame_width=500, frame_height=400))

mpd.show_versions()

MovingPandas 0.16.0

SYSTEM INFO
-----------
python     : 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:16:53) [MSC v.1929 64 bit (AMD64)]
executable : H:\miniconda3\envs\mpd-ex\python.exe
machine    : Windows-10-10.0.19045-SP0

GEOS, GDAL, PROJ INFO
---------------------
GEOS       : None
GEOS lib   : None
GDAL       : 3.5.0
GDAL data dir: None
PROJ       : 9.0.0
PROJ data dir: H:\miniconda3\pkgs\proj-9.0.0-h1cfcee9_1\Library\share\proj

PYTHON DEPENDENCIES
-------------------
geopandas  : 0.13.0
pandas     : 2.0.1
fiona      : 1.8.21
numpy      : 1.23.5
shapely    : 1.8.2
rtree      : 1.0.1
pyproj     : 3.3.1
matplotlib : 3.7.1
mapclassify: 2.4.3
geopy      : 2.3.0
holoviews  : 1.14.9
hvplot     : 0.8.3
geoviews   : 1.9.6
stonesoup  : 0.1b12
In [2]:
gdf = read_file('../data/geolife_small.gpkg')
traj_collection = mpd.TrajectoryCollection(gdf, 'trajectory_id', t='t')

Clipping a Trajectory¶

In [3]:
xmin, xmax, ymin, ymax = 116.365035,116.3702945,39.904675,39.907728
polygon = Polygon([(xmin,ymin), (xmin,ymax), (xmax,ymax), (xmax,ymin), (xmin,ymin)])
polygon_gdf = GeoDataFrame(pd.DataFrame([{'geometry':polygon, 'id':1}]), crs=31256)

my_traj = traj_collection.trajectories[2]
intersections = my_traj.clip(polygon)
    
print("Found {} intersections".format(len(intersections)))

Found 1 intersections
In [4]:
ax = my_traj.plot()
polygon_gdf.plot(ax=ax, color='lightgray')
intersections.plot(ax=ax, color='red', linewidth=5)
Out[4]:
<Axes: >

Clipping a TrajectoryCollection¶

Alternatively, using TrajectoryCollection:

In [5]:
clipped = traj_collection.clip(polygon)
clipped
Out[5]:
TrajectoryCollection with 2 trajectories
In [6]:
clipped.plot()
Out[6]:
<Axes: >

Computing intersections¶

In [7]:
polygon_feature = {
    "geometry": polygon,
    "properties": {'field1': 'abc'}
}
In [8]:
my_traj = traj_collection.trajectories[2]
intersections = my_traj.intersection(polygon_feature)
intersections
Out[8]:
TrajectoryCollection with 1 trajectories
In [9]:
intersections.to_point_gdf()
Out[9]:
id sequence trajectory_id tracker geometry intersecting_field1
t
2009-02-04 10:43:04.222205 3225 773 3 2 POINT (116.36799 39.90467) abc
2009-02-04 10:43:05.000000 3226 774 3 2 POINT (116.36798 39.90471) abc
2009-02-04 10:43:07.000000 3227 775 3 2 POINT (116.36795 39.90480) abc
2009-02-04 10:43:09.000000 3228 776 3 2 POINT (116.36793 39.90487) abc
2009-02-04 10:43:11.000000 3229 777 3 2 POINT (116.36794 39.90495) abc
... ... ... ... ... ... ...
2009-02-04 10:48:14.000000 3390 938 3 2 POINT (116.36535 39.90589) abc
2009-02-04 10:48:15.000000 3391 939 3 2 POINT (116.36525 39.90589) abc
2009-02-04 10:48:16.000000 3392 940 3 2 POINT (116.36516 39.90589) abc
2009-02-04 10:48:17.000000 3393 941 3 2 POINT (116.36507 39.90589) abc
2009-02-04 10:48:17.399995 3393 941 3 2 POINT (116.36504 39.90589) abc

170 rows × 6 columns

In [ ]: